#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use db_parameters;
use ortho_species;
use die_error;
use run_cmd;

use DBI;


print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_kaks_align

	Description:
		Enter Kaks and alignment data for pairwise relationships
			within clades of ortholog_sets
		Includes both orthologues and Orphan inparalogs (with pseudos)
		
	Populates DB tables:
		orthologs.kaks_align
	Relies on:
		ens_id
		ortholog_sets
		(especially after phyop orthologs resolved)

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	get_prot_id_to_gene_id_pi_map
#
#
sub get_prot_id_to_gene_id_pi_map($\%\%)
{
	print STDERR "\tGet prot_id to gene_id and PI mapping\n";
	my ($dbh, $prot_id_to_gene_id, $prot_id_to_pi) = @_;
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT 
            prot_id, 
            gene_id,
            pi
        FROM
            ens_id
        WHERE
            species = '$ortho_name1' or
            species = '$ortho_name2';
PL/SQLCMD
	my  $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd."\n";
	die unless scalar @$data;
	for my $line(@$data)
	{
		$prot_id_to_gene_id->{$line->[0]} = $line->[1];
		$prot_id_to_pi->{$line->[0]} = $line->[2];
	}
}


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	get_ortholog_prot_id_pairs
#
#
sub get_ortholog_prot_id_pairs($$)
{
	my ($dbh, $protocol_id) = @_;
	print STDERR "\tGet intra clade protein identifier pairs\n";
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT 
            e1.prot_id AS prot_id1,
            e2.prot_id AS prot_id2
        INTO 
            TEMP t_intra_clade_pairs 
		FROM 
            ens_id e1 JOIN 
			orthologs.ortholog_sets os1 using (gene_id) JOIN 
			orthologs.ortholog_sets os2 USING (protocol_id, ortholog_id) JOIN
            ens_id e2 on (os2.gene_id = e2.gene_id)
        WHERE 
            protocol_id = $protocol_id AND
            e1.prot_id < e2.prot_id;
PL/SQLCMD
	$dbh->do($sql_cmd) or die die $dbh->errstr.".";

	
	print STDERR "\tGet closest same species protein identifier pairs\n";
$sql_cmd =<<"PL/SQLCMD";
   INSERT INTO t_intra_clade_pairs
      SELECT 
            e1.prot_id AS prot_id1,
            e2.prot_id AS prot_id2
         FROM 
            ens_id e1 JOIN 
            orthologs.closest_same_species css using (gene_id) JOIN 
            ens_id e2 on (css.same_species_gene_id = e2.gene_id)
         WHERE 
            protocol_id = $protocol_id AND
            e1.prot_id < e2.prot_id;
PL/SQLCMD
	$dbh->do($sql_cmd) or die die $dbh->errstr.".";


	print STDERR "\tGet closest other species protein identifier pairs\n";
$sql_cmd =<<"PL/SQLCMD";
   INSERT INTO t_intra_clade_pairs
      SELECT 
            e1.prot_id AS prot_id1,
            e2.prot_id AS prot_id2
         FROM 
            ens_id e1 JOIN 
            orthologs.closest_other_species cos using (gene_id) JOIN 
            ens_id e2 on (cos.other_species_gene_id = e2.gene_id)
         WHERE 
            protocol_id = $protocol_id AND
            e1.prot_id < e2.prot_id;
PL/SQLCMD
	$dbh->do($sql_cmd) or die die $dbh->errstr.".";


	print STDERR "\tRetrieve protein identifier pairs\n";
$sql_cmd =<<"PL/SQLCMD";
   SELECT
         DISTINCT ON (prot_id1, prot_id2)
         prot_id1, prot_id2
      FROM 
         t_intra_clade_pairs;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd);
	$dbh->do("DROP TABLE t_intra_clade_pairs;");
	return $data or die_error DBI::errstr.$sql_cmd."\n";
}


#
# connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
$dbh->{RaiseError} = 0;




#
#	Get prot_id to gene_id mapping
#
my (%prot_id_to_gene_id, %prot_id_to_pi);
get_prot_id_to_gene_id_pi_map($dbh, %prot_id_to_gene_id, %prot_id_to_pi);


for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;

	# hack for ensembl / inparanoid
	#
	#	next unless $protocol_id == $curr_protocol_id;
	print STDERR "\tOrthology protocol $protocol_name...\n";
	$dbh->do("DELETE FROM orthologs.align_kaks WHERE protocol_id = $protocol_id;");


	my $species_pair	= $ortho_name1.'_'.$ortho_name2;
	my $file_ids		= "$dir_pipeline_temp/prot_id.pairs";
	my $file_kaks		= "$dir_orthologs_output/dnds_res.both_filtered";

	# hack for ensembl / inparanoid
	#
	# if ($protocol_name eq 'inparanoid' or $protocol_name eq 'ensembl')
	# {
	# 	$file_kaks		= "$dir_orthologs_output/merged.kaks_res" ;
	# }

	
	my $file_kaks_temp	= "$dir_pipeline_temp/orthologs.kaks";
	print STDERR "\tIdentifiers at $file_ids      \n";
	print STDERR "\tKaKs        at $file_kaks     \n";
	print STDERR "\tKaks temp   at $file_kaks_temp\n";

	
	#
	#	get prot_id pairs
	#
	print STDERR "\tGet protein identifier pairs...\n";
	my $prot_id_pairs = get_ortholog_prot_id_pairs($dbh, $protocol_id);
	{
		print "\t", scalar @$prot_id_pairs, " prot_id pairs\n";
		open_or_die (*IDS, ">$file_ids");
		print IDS (join ("\t",  @$_). "\n") for (@$prot_id_pairs);
	}

	#
	#	get kaks values
	#
	print STDERR "\n\tGet matching entries from [$file_kaks]...\n";
	run_cmd (<<"bourneshell_CMD");
	$dir_lg_share_bin/seq_pairs_matching
		--file1in $file_kaks
		--file2in $file_ids
		--shared_in_full $file_kaks_temp
		--counts $dir_pipeline_ds_temp/summary/orthologs.pairwise.counts
		--verbose
bourneshell_CMD
	

	open_or_die(*KAKS, $file_kaks_temp);

	
	#
	#	Enter into panda
	#
	print STDERR "\t\tEntering into panda...\n";
	$dbh->do("COPY orthologs.align_kaks(".
					"gene_id1, gene_id2, prot_id1, prot_id2, pi1, pi2, dn, ds, dnds, count_n, ".
					"count_s, dn_stderr, ds_stderr, kappa, lnl, tau, error_str, align_start, ".
					"align_len, aligned_positions, coverage, e_value, len1, len2, ".
					"aligned_residues, percent_id, bitscore, gaps1, gaps2, protocol_id".
					") FROM STDIN") or die DBI::errstr;

	# initialize dots
	my $dots_total = 400000;
	my $dots_expect_lines = 4;
	my $dots_counter = 0;
	my $dots_every = $dots_total / ($dots_expect_lines * 50) + 1;
	print STDERR "\t[";
	
	while (<KAKS>)
	{
		# Update dots
		if (++$dots_counter % $dots_every == 0) 
		{
			my $dots_line = int($dots_counter / $dots_every);
			print STDERR "\n\t" if ($dots_line && $dots_line % 50 == 0); 
			print STDERR ".";
		}
		
		
		chomp;
		my ($prot_id1,  $prot_id2, @data)  = split (/\t/, $_, -1);

		$dbh->pg_putline(join ("\t",  
								$prot_id_to_gene_id{$prot_id1},
								$prot_id_to_gene_id{$prot_id2},
								$prot_id1,
								$prot_id2,
								$prot_id_to_pi{$prot_id1},
								$prot_id_to_pi{$prot_id2},
								@data,
								$protocol_id)."\n");
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy  or die DBI::errstr."\nEND COPY\n";

	# Finish dots
	print STDERR "]\n";
	
}
$dbh->do("VACUUM  orthologs.align_kaks;");


print STDERR "\tCompleted\n";


log_pipeline_stage($dbh, $curr_protocol_id, 39, 'table kaks_align', $start_time);
$dbh->disconnect();

